import delimited "${tmp}/labexkmean.csv", delim(",") encoding(utf8) clear
	keep labex kmeans
	gduplicates drop
	
	rename labex mainlabex
	merge 1:m mainlabex using  "${tmp}/Regfile_nafZE", nogen keep(2 3)
	qui su kmean
	replace kmean = `r(max)' + 1 if mi(kmean)
	

	xtset gp year
	drop _m
	
	merge m:1 ZE2010 using "${data}/Utils/ZE_regIDF_mapping"
	

	// define sample 
	egen mlabexid = group(mainlabex)
	qui su mlabexid
	replace mlabexid = r(max) + ZE2010 if mi(mlabexid)
	

	reghdfe ln_eng_sbrut E_post11  if $cond_accept, a(ZE2010#year  gp)  cl(mlabexid)	
	bys gp: egen insampeng = max(e(sample))
	keep if insampeng

	
	local row0 = "& Static Estimate & Pre Trends test & Obs. \\"

	local row1 = "Sample restrictions \\"
	preserve
	drop if IDF == 1
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year") clust(mlabexid) 
	local row2 = "1. Removing Paris & `=string(e(coef), "%9.4f")'`=e(star)' & `=e(pt)' & `=e(Nobs)' obs \\"
	local row3 = "& (`=string(e(se), "%9.4f")') & & (`=e(Npairs)' pairs) \\ "
	restore
	
	preserve
	drop if naf2d == 72
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year") clust(mlabexid) 
	local row4 = "2. Removing R&D sector & `=string(e(coef), "%9.4f")'`=e(star)' & `=e(pt)' & `=e(Nobs)' obs \\"
	local row5 = "& (`=string(e(se), "%9.4f")') & & (`=e(Npairs)' pairs) \\ "
	restore
	
	preserve
	drop if naf2d == 20
	drop if naf2d == 21
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year") clust(mlabexid) 
	local row6 = "3. Removing chemistry and pharma sectors & `=string(e(coef), "%9.4f")'`=e(star)' & `=e(pt)' & `=e(Nobs)' obs \\"
	local row7 = "& (`=string(e(se), "%9.4f")') & & (`=e(Npairs)' pairs) \\ "
	restore	
		
	xtset gp year
	xtreg D.ln_eng_sbrut if year <= 2010, fe
	predict temp, xbu
	bys gp: egen gp_croiss = max(temp)
	gen gp_trend = (year - ${NORMYEAR})*gp_croiss
	gen gp_trend_sq = gp_trend^2
	gen gp_trend_cub = gp_trend^3
	gen gp_trend_quad = gp_trend^4

	local row8 = "Controlling for trends"
	preserve
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year") controls(gp_trend) clust(mlabexid) 
	local row9 = "1. Pre-reform dep. var. linear trend & `=string(e(coef), "%9.4f")'`=e(star)' & `=e(pt)' & `=e(Nobs)' obs \\"
	local row10 = "& (`=string(e(se), "%9.4f")') & & (`=e(Npairs)' pairs) \\ "
	restore	
	
	preserve
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year") controls(gp_trend gp_trend_sq gp_trend_cub gp_trend_quad) clust(mlabexid) 
	local row11 = "2. Pre-reform dep. var. trend 4th deg. poly & `=string(e(coef), "%9.4f")'`=e(star)' & `=e(pt)' & `=e(Nobs)' obs \\"
	local row12 = "& (`=string(e(se), "%9.4f")') & & (`=e(Npairs)' pairs) \\ "
	restore	
	

	gen gp_trendt = gp_trend*log(1+exposure)
	gen gp_trendt2 = gp_trend_sq*log(1+exposure)
	gen gp_trendt3 = gp_trend_cub*log(1+exposure)
	gen gp_trendt4 = gp_trend_quad*log(1+exposure)
	gen lexposure = log(1+exposure)
	
	preserve
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year") controls(gp_trend*) clust(mlabexid) 
	local row13 = "3. Pre-reform dep. var. trend 4th deg. poly x shock exposure & `=string(e(coef), "%9.4f")'`=e(star)' & `=e(pt)' & `=e(Nobs)' obs \\"
	local row14 = "& (`=string(e(se), "%9.4f")') & & (`=e(Npairs)' pairs) \\ "
	restore	
	
	forvalues i = 1/12 {
		gen km`i' = kmeans == `i'
		gen kmI`i' = km`i'*year
	}
	
	preserve
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year") controls(kmI* km*) clust(mlabexid) 
	local row15 = "4.Flexible discipline trend & `=string(e(coef), "%9.4f")'`=e(star)' & `=e(pt)' & `=e(Nobs)' obs \\"
	local row16 = "& (`=string(e(se), "%9.4f")') & & (`=e(Npairs)' pairs) \\ "
	restore		

	
	gen pt = prox*year
	preserve
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year") controls(pt) clust(mlabexid) 
	local row17 = "5.Flexible science prox. trend & `=string(e(coef), "%9.4f")'`=e(star)' & `=e(pt)' & `=e(Nobs)' obs \\"
	local row18 = "& (`=string(e(se), "%9.4f")') & & (`=e(Npairs)' pairs) \\ "
	restore		
	
	local row19 = "Alternative inference (p-value) \\"
	
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year")  clust(mlabexid) 
	local row20 = "1. Baseline & `=e(pm)' & `=e(pt)' \\"
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year")  clust(ZE2010) 
	local row21 = "2. Clustered at the CZ level& `=e(pm)' & `=e(pt)' \\"	
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year")  clust(gp) 
	local row22 = "3. Clustered at the unit level &`=e(pm)' & `=e(pt)' \\"	
	qui prog_gentable E_2005, estimator("reghdfe")  dep(ln_eng_sbrut) fe("gp ape_n#year ZE2010#year")  clust(ape_n ZE2010) 
	local row23 = "4. Twoway clustering & `=e(pm)' & `=e(pt)' \\"		
	

	
	
use "${tmp}/Proximity_naflabex", clear

	merge m:1 labexid using "${data}/Utils/labexcom", nogen keep(1 3)

	replace depcom = 75056 if floor(depcom/1000) == 75
	replace depcom = 13055 if inrange(depcom,13201,13216)
	replace depcom = 69123 if inrange(depcom,69381,69389)

	merge m:1 depcom using "${data}/Utils/depcom_ZE_mapping", nogen keep(1 3)

	keep labexid ZE2010
	duplicates drop
	drop if ZE == ""
save "${tmp}/ZE_labex", replace



forvalues r = 1/10000 {
 qui {
 use labexid dotation accept using "${tmp}/main_labex_with_ctrf", clear
	merge 1:1 labexid using "${tmp}/ZE_labex", keep(3) nogen
	
	preserve
		keep dotation 
		drop if dotation == 0  | dotation == .
		gen rdn = _n
		save "${tmp}/rdn_dotation", replace
	restore
	egen nb_rejected = sum(accept == 0)
	
	keep ZE2010 accept dotation labexid nb_rejected
	gen foo = runiform()
	sort foo
	su nb_rejected
	sort foo
	gen reject = _n <= `r(mean)'
	sort reject foo
	
	gen rdn = _n if reject == 0
	keep ZE2010 accept dotation labexid rdn
	rename dotation dotation_ori
	merge m:1 rdn using "${tmp}/rdn_dotation"
	rename dotation dotation_perm
	recode dotation* (missing = 0)
	
	
	keep  labexid ZE2010 dotation*
save "${tmp}/labexperm", replace
	
use "${tmp}/matrix_issn_labex", clear
	joinby issn using "${tmp}/matrix_ape_issn"
	bys ape labex: egen share = sum(share_labex_issn * share_ape_issn)
	keep labex share ape
	duplicates drop
		
	rename share prox

	recode prox (missing = 0)
		
	merge m:1 labexid using "${tmp}/labexperm", nogen keep(3)
		
	gen exposure_perm = prox*dotation_perm
	gcollapse (sum) exposure_perm , by(ZE2010 ape) fast
		
	destring ZE, replace force
save "${tmp}/Proximity_naflabex_perm", replace

use "${tmp}/Regfile_nafZE", clear
	merge m:1 ZE2010 ape using "${tmp}/Proximity_naflabex_perm", keep(1 3) nogen
		
	replace exposure_perm = 0 if exposure_perm == .
	forvalues year = 2004/2019 {
		gen  EPerm_`year' = ln(1+exposure_perm) * (year == `year')
	}
	gen EPerm_post11 = ln(1+exposure_perm) * (year >= 2011)
	keep if inrange(year, 2004, 2019)
	drop EPerm_${NORMYEAR} 

	reghdfe ln_eng_sbrut   E_post11  if $cond_accept, a(ZE2010#year  gp)  cl(gp)	
	bys gp: egen insampeng = max(e(sample))
	
	reghdfe ln_eng_sbrut   EPerm_post11 if $cond_accept & insampeng, a(ZE2010#year ape_n#year gp) 	
	local coef_stat_fe = _b[EPerm_post11]
	reghdfe ln_eng_sbrut   EPerm_???? if $cond_accept & insampeng, a(ZE2010#year ape_n#year gp) 
	lincom 1/6*(EPerm_2004+EPerm_2005+EPerm_2006+EPerm_2007+EPerm_2008+EPerm_2009) 
	local pre_trend_fe = `r(estimate)'	
	
	clear
	gen coef_stat_fe = .
	gen iter = .
	gen coef_pretrend_fe = .
	set obs 1
	replace iter = `r'
	replace coef_stat_fe = `coef_stat_fe'
	
	replace coef_pretrend_fe = `pre_trend_fe'
	
	if `r' != 1 {
		append using "${tmp}/perm_result"
	}
	save "${tmp}/perm_result", replace
}
if mod(`r',1) == 0 {
	disp in red "Iteration number: `r'"
	}
}
		
	
use "${tmp}/Regfile_nafZE", clear
		
		keep if inrange(year, 2004, 2019)
		drop if  inrange(naf2d,64,64)  |  naf2d == 70  /* | naf2d == 77 */

		reghdfe ln_eng_sbrut E_post11  if $cond_accept  , a(ZE2010#year ape_n#year gp) 
		local coef_stat_fe = _b[E_post11]
		reghdfe ln_eng_sbrut    E_???? if $cond_accept , a(ZE2010#year ape_n#year gp) 
		lincom 1/6*(E_2004+E_2005+E_2006+E_2007+E_2008+E_2009) 
		local pre_trend_fe = `r(estimate)'		
clear
	gen coef_stat_fe = .
	gen coef_stat_nofe = .
	gen iter = .
	gen coef_pretrend_fe = .
	set obs 1
	replace iter = 0
	replace coef_stat_fe = `coef_stat_fe'
	replace coef_pretrend_fe= `pre_trend_fe'


	forvalues i = 2004/2019 {
		if `i' != ${NORMYEAR} {
			gen coef_`i' = _b[E_`i']
		}
	}
	append using "${tmp}/perm_result"

	
		su coef_stat_fe if iter == 0
		local mfe = r(mean)
		count if abs(coef_stat_fe) > abs(`mfe') 
		local c1 = r(N)
		su coef_pretrend_fe if iter == 0		
		local ptfe = r(mean)
		count if abs(coef_pretrend_fe) > abs(`ptfe')
		local row24 = "5. Randomization based inference & [p=0.`=round(`c1'/10000, 0.001)'] & [p=0.`=round(`r(N)'/10000, 0.001)'] \\"		
		
			
			file open myfile using "${outpath}/tableA6.tex", write text replace
	forvalues i = 0/24 {
		file write myfile " `row`i''" _n
	}
	file close myfile
	
		/*
 		gen foo = (coef_stat_fe - `mfe')^2
		egen stderr = sum(foo/10000)
		su stderr
		twoway  (hist coef_stat_fe if iter != 0, col(red%50) xline(`mfe', lcol(red)))

	